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We study numerically the growth of a crack in an elastic medium under the influence 
of a travelling Shockwave. We describe the implementation of a fast algorithm which 
is perfectly suited for a data parallel computer. Using large scale simulations on the 
Connection Machine we generate cracks with more than 10000 sites on a 1024x1024 
lattice. We show that the resulting patterns are fractal with a fractal dimension that 
depends on the chosen breaking criterion and varies between 1. and 2. 



1. Introduction 

How does a solid body break under an externally applied load? This is an important 
question for researchers in many fields and has strong implications on e.g. mate- 
rials science, engineering or geophysics. It has been studied for a long time quite 
extensively using numerical, experimental and analytical methods.^ 

If one considers the solid to be a linear elastic medium the formulation of the 
initial question can be made in terms of the Lame equation 

(l-2i/)Aw + V(V-w) = (1) 

which describes the displacement u of a small volume element in an elastic material 
from its equilibrium position. The Poisson ratio u i&a, material dependent parameter 
which has the following meaning: If one applies a uniaxial force to an elastic bar 
of length L and cross section W x VF, it will not only change its length by Ai in 
the direction of the force, but it will also change its width by /S.W in the direction 
perpendicular to the force. The Poisson ratio then defines the relative amount of 
change 



AL/L 



(2) 



Due to thermodynamical reasons v is bounded between — 1 < < 0.5. Concrete 
for instance has a Poisson ratio v « 0.2. 
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A crack in such a system can be described as an additional force free surface 
which thus obeys the boundary condition 

g-n = (3) 

where a is the stress tensor and n is a surface vector. The crack grows if the stress 
parallel to the crack surface (T|| is larger than a certain material dependent critical 
stress ac- Since no first principle law for the normal growth velocity is known one 
assumes the general behavior 

Vn oc (a|| - ac)'' (4) 

where r] is often simply set to unity. The equations (1), (3) and (4) formulate the 

growth of a crack as a moving boundary problem which is far more difficult to solve 
than the Lame equation itself. Yet, there is one important property of real material 
missing, which is "disorder". Microscopically disorder means deviations from the 
perfect crystal structure of the clastic material due to vacancies, dislocations or grain 
boundaries. But macroscopically these imperfections are simply reduced to spatial 
randomness of the material properties like Poisson ratio and breaking threshold. 

In numerical simulations for crack growth one often discretizes the elastic 
medium on a lattice, for instance with a finite clement scheme.^ According to the 
chosen boundary conditions and the externally applied load one relaxes the system 
to equilibrium. Then one picks one or several bonds according to a given rule - 
for instance the bond with the highest load and breaks it. This defines a new 
boundary and therefore one has to resolve the whole problem again. This procedure 
is then repeated until a crack of desired size is grown. The disorder is often put into 
the simulations by considering lattice parameters (like bond strength and breaking 
threshold) that vary from site to site. Since the system has enough time to relax 
to full equilibrium before the crack can grow, such a procedure is only capable of 
describing quasi static processes. 

On the other hand there are phenomena — like explosions, shattering of glass 
or shock waves - which should not be treated in a static approximation and in 
which the system is not able to relax to equilibrium before breaking a bond, but 
in which the time to relax the system is comparable to or larger than the time one 
needs to propagate the crack. Such a situation usually results in cracks with many 
sidebranches growing behind the shock front since the internal energy cannot be 
dissipated fast enough.^ In the following we are going to study numerically such a 
process in which one obviously has to take into account the dynamical behavior of 
the elastic medium. 

The organization of this paper is thus as follows: in Sect. 2 we introduce the 
general model, whose implementation on the Connection Machine is described in 
Sect. 3. In Sect. 4 we describe the details of the simulation and in Sect. 5 we present 
some results. 

2. The Model 

As a model for the elastic material we consider a triangular network of Hookean 
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springs connecting points of mass m. A triangular network is necessary for the 
simulation since a simple quadratic lattice does not have any shear modulus and 
thus can be deformed arbitrarily under shear load. This model is known as "central 
force model" ^ since the Hookean springs are isotropic. The boundary sites of this 
network are kept fixed in space. On each site there are two continuous degrees of 
freedom, which are the coordinates of the displacement Ux and Uy of this site from 
its equilibrium position rg. Since we want to study the dynamical behavior of this 
network we have to determine the time evolution of the displacements which is gov- 
erned on each lattice site by Newton's equation. Following a suggestion of Chopard^ 
we use a discrete time Hamilton formalism to express the dynamic behavior. The 
Hamiltonian of this system is given by 



h{p^.. .p^,r, . . . r^) = ^ + - ^ (r„r,.) (5) 



where p. and are the momentum and position of site i, and Uij is the interaction 

energie between two lattice sites i and j. Uij is nonzero only between nearest 
neighbor sites and since we use Hookean springs it is simply a harmonic potential 

Uij = ilu -rj\- = % {\ui +Loi- Uj -Loj\- a)^ 

ki 2 
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where kij is the coupling constant between the sites i and j. dr^ = Tq^ —Loj is 
the vector between their equilibrium positions and a is the equilibrium length of 
the connecting spring (please note, that dr^j does not mean differentials). In our 
simulation we arc going to set a = while we keep dr^ = 1. This corresponds to 
the case of a prestretched network and can be compared for instance with the skin 
of a drum. The discretized versions of Hamilton's equations are 

p_^[t)-p_^{t-l) = -— 

where we use vector derivatives as symbols for the corresponding gradients. The 
time evolution of the displacement field can finally be written as 



1 dH 

W,(t+l)-2u,(i)+U,(t-l): 



rrii du^(t) 



(8) 



j=NN{i) 



where S^j = Ui{t) —Uj{t) + dr^j . This equation (8) defines the updating rule, which 
relates the displacements at time t + 1 to the displacements at times t and t — 1. 
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Although the chosen dynamics seems to be rather crude as compared to molecular 
dynamics, one can actually show that it conserves the total momentum and some 
kind of total "energy" . Moreover we show that it conserves the essential features 
we need for the generation of cracks: A wave packet that is imposed on the lattice 
will travel with only minor changes in shape with a definite velocity through the 
lattice. 

3. The Growth Model using Fortran 90 

Since the same updating rule (8) is applied to all lattice sites at each time step and 
on the other hand the topology of the underlying lattice is, in contrast to molecular 
dynamics, not changed dnc to rearrangements of particles single instruction multi- 
ple data (SIMD) machine like the CM is the appropriate computer architecture for 
this problem: Each lattice site is mapped onto one virtual processor and all sites 
are updated in parallel. The only inter processor communication is required for the 
calculation of the right hand side of eqn. (8). But because the triangular lattice 
structure can be mapped onto a square lattice with next nearest neighbor interac- 
tions into one direction, a nearest neighbor grid communication using CSHIFTS is 
sufficient and no general router communication is required, which makes this lattice 
model fast and efficient. 

In the actual implementation of the program we chose the following data layout 
for the main variables: Each processor has to store the displacements at time t and 
< — f , UT and UTMl, which are represented as complex numbers. Since we intend 
to simulate a "disordered" system, we assign to each spring a different, randomly 
chosen coupling constant. Therefore we keep on each lattice site the six couplings 
K to all neighbors. By doing this we waste some memory space because each kij 
is stored twice - on site i and on site j - but we save computer time by avoiding 
unnecessary communication. The data layout for the main variables is thus as 
follows 

COMPLEX , ARRAY (NXY, NXY) : : UT, UTMl 

REAL , ARRAY (Z, NXY, NXY) :: K 
CMF$ LAYOUT UT ( : NEWS , : NEWS) 
CMF$ LAYOUT UTMl ( : NEWS , : NEWS) 
CMF$ LAYOUT K ( : SERIAL , : NEWS , : NEWS ) 

Listing 1. 



The LAYOUT directive for K is used to group the couplings of one site as a serial di- 
mension onto one processor. It is necessary, because otherwise the compiler would 
spread the couplings over all virtual processors which results in unnecessary com- 
munication for the force calculation. Using this data layout a single step of the 
updating rule (8) is programmed in a straightforward manner in CM Fortran 

C GET VECTORS TO NEIGHBORS 

C COMPUTE NEW PARTICLE POSITION AND MOMENTUM 

C NN is a help field to store r_i-r_j . It has the same layout as 
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C ut and utml. 



NN(1,:,:) = CSHIFT (UT , 1 , 1) 

NN(2,:,:) = CSHIFT(UT,2 , 1) 

NN(4,:,:) = CSHIFT (UT, 1,-1) !get displacements 

NN(5,:,:) = CSHIFT (UT, 2,-1) !of nearest neighbors 

NN(3,:,:) = CSHIFT(NN(4, : , : ) ,2, 1) 

NN(6,:,:) = CSHIFT(m(l , : , : ) ,2, -1) 

NW(1,:,:) = WN(1, : , :)-UT+DR_l 

NN(2,:,:) = NN(2, : , : )-UT+DR_2 

NN(3,:,:) = NN(3, : , : ) -UT+DR_3 

MN(4,:,:) = NN (4 , : , : ) -UT+DR_4 {calculate difference vector 

NN(5,:,:) = NN(5, : , : )-UT+DR_5 
NN(6,:,:) = NN(6, : , : )-UT+DR_6 



UTMl = SUM(NN*(K*(1.-A1/ABS(NN))) ,DIM=1)/M + 2.*UT - UTMl 

Listing 2. 



The CSHIFT commands are used to communicate the displacements between "neigh- 
boring" processors. By carefully reusing already shifted data it is of course possible 
to get the displacements from the six nearest neighbors on the triangular lattice with 
only six CSHIFTS. Since we are not using full next nearest neighbor communication 
the use of stencil operations does not seem useful at this point. The global SUM along 
the first dimension computes the total force exerted onto each site by its neighbors. 
Because the first dimensions of the coupling constant array k and the vectors to the 
nearest neighbors NN are laid out onto the same processor as a : serial dimension, 
the computation of the SUM does not require any communication. 

Unfortunately this simple formulation docs not give optimal performance. The 
compiler allocates and deallocates unnecessary temporary fields and even generates 
general CM_send router communication, which in fact uses 43% of the total CPU 
time! To obtain a much bettor performance we coded the code fragment shown 
in Listing 2 completely in PARIS (PARallel Instruction Set). This allows to fully 
control the memory allocation, the communication and to make efficient use of 
pipelined commands like CM_f _sub_constjnult_always. Now, most of the CPU 
time ,38% - is used for the calculation of the distance between neighboring lattice 
sites |r — r^ l which is coded as CM_f _c_abs_2. Now, the NEWS communication part 
is negligible and sums up to 5.8% of the total CPU time. 

These improvements result in a speedup of a big factor of five as compared to 
the straightforward implementation. One obtains on 8K processors of the previously 
described CM2 an update rate of 1.1 miUions of updates per second (MUPS) and 
a speed of 110 MFlops. As a comparison, one can obtain with typical spin cellular 
automata using multispin coding techniques more than 1000 MUPS on one processor 
of a NEC-SX35 and on a CM2-16K a Q2R cellular automaton runs at 1600 MUPS.*^ 
For a full molecular dynamics simulations on an CM200-8K Hedman and Laaksonen 
achieve about 0.2 MUPS'^ and for MD simulations on a CM2-16K Mel'cuk et. al. 
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obtained an update rate of 4.5 KUPS.^ 

4. The simulation 

All simulations arc performed on a 1024 x 1024 lattice and we use between 3000 
and 5000 timesteps. At the beginning of the simulation the couplings are chosen at 
random out of a uniform distribution with mean value ko = 0.005 and a width of 
typically 50%. To initiate the crack growth we break all bonds between site rg and 
its neighbors. Afterwards an initial pulse is imposed on the center of the lattice: 
If Tq is the central site in the lattice, then its nearest neighbors . . .r^ are dis- 
placed radially outward for 100 lattice units while keeping all other sites fixed. This 
displacement seems very large at first sight, but since we use a harmonic potential 
between the sites the actual size of the initial displacement is not relevant. We 
only wanted to make sure that this perturbation is much larger than the "thermal" 
motion unduced by chosing random coupling constants. At time t = all sites are 
released and the system is free to evolve. After every other time step one looks 
for the bond with the largest elongation Im- Then one determines all bonds whose 
elongations I are larger than a ■ Im - where a is an adjustable parameter - and 
which lie on the surface of the already existing crack. All those bonds are broken 
by setting the corresponding coupling constant kij to zero. Here one has to notice 
that this breaking rule does not require any communication and can be performed 
completely in parallel since both sites i and j which are connected by such a bond 
will clear their own copy of kij . 

Thus wc consider a relative breaking threshold rather than an absolute one, 
which has the following reason; In an absolute breaking threshold one would break 
all bonds whose elongation is larger than some fixed critical length Ic- On the other 
hand the amplitude of the outgoing wave packet is decreasing with the distance from 
the center. So, when the wave packet has initially an amplitude that is larger than 
the threshold, the outgoing wave will break all bonds it reaches until its amplitude 
has dropped below the threshold and from then on no further bond will be broken. 
Thus, one only creates a structureless isotropic hole in the center of the system. 

A somewhat similar model has been studied by Louis et. al.^ They try to solve 
a quasi static problem, but perform only a few relaxation steps to find the equilib- 
rium state of the network. However, since they use an overrelaxation scheme the 
relaxation of their system towards equilibrium has another dynamical interpretaion 
than the iterative method I use. Another similarity is given in the breaking rules. 
Louis et. al. pick bonds that are to be broken with a probability that is proportional 
to the bond length oc |rj — Lj\- 

5. Results 

To demonstrate that the dynamics (8) produces reasonable results we first consider 
the case of a smooth wave packet travelling through a system without disorder and 
without breaking. Therefore we applied not a singular pulse to the network but 
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rather a smooth wave packet. The initial radial displacements of the central sites 
are 

u{t) =Uq-(^- cos ^27r • ^ ^ (9) 

for < t < r and u{t) = for < > r while all other sites are free to move. The 
period is chosen to be r = 100. 




-30.00 I • ^ • ^ • ^ • ^ • ^ • ^ • ^ • 1 

-64.00 -48.00 -32.00 -16.00 0.00 16.00 32.00 48.00 64.00 

X coordinate 



Fig. 1. Cut through a Shockwave after 200 timesteps. 

Fig. (1) shows a cut along the x axis through this wave packet at time t = 200. One 

indeed recovers the original pulse plus a further minimum which is due to the fact 
that we keep the displacements of the central sites for t > r fixed at zero. Because 
of their inertia the neighboring sites keep vibrating which results in the second 
minimum. If one measures the velocity of the maximum of the wave packet one 
obtains a value which agrees with the analytical expression for the group velocity 
of a circular wave with frequency w = 27r/r on a triangular lattice. 

After having confirmed that the travelling wave shows reasonable behavior we 
restrict ourselves again to the case of a singular pulse, which corresponds to the 
case T = 2. 

Fig. 2. Crack patterns generated for (from left to right) a) a = 1, b) a = 0.98 and c) a = 0.95. 

In figs (2.a-2.c) we show examples of cracks which were generated for different 
breaking thresholds: a = 1 (10934 sites) - which means that only the bond with 
the largest elongation is broken - a = 0.98 (16038 sites) and a = 0.95 (35837 sites). 
The number in brackets are the number of sites that are "connected" by broken 
bonds. All cracks show a starlike and fractal structure. For decreasing a they 
become more and more ramified. This is easily understood since with decreasing 
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a many more bonds are eligible to be broken. The delta-like excitation we impose 
initially on the lattice is highly non-periodic and therefore produces waves with 
many different frequencies, althoug lower amplitude. For small enough a. many of 
these waves can contribute to the growth of the crack. But one also has to take into 
account another effect. Since the lattice is prestretched each bond that is broken is a 
source for another spherical wave travelling away from this point. For small enough 
a also those waves can break bonds and therefore can lead to an avalanche-type 
growth of the crack. Thus, one can distinguish two different regimes: For a close 
to unity the crack grows mainly at the tips at the outer branches which coincide 
with the front of the shockwavc. The sidcbranchcs behind the shockfront remain 
inactive. For smaller a also the tips behind the shockfront continue to grow and 
split and eventually the crack becomes space filling. 

Another fact to be noticed is that for decreasing a the lattice structure becomes 
more and more dominant and eventually the cracks grow into a structure with sixfold 
symmetry. 

To be more quantitative, we study the dependence of the number of broken 

bonds N on the radius of the cluster, which is measured in terms of the radius of 
gyration Rq. Rg describes the average distance of all broken bonds from the center 
of mass T.Qj^ of the crack 

N 

Y.(^r-^-CMf- (10) 





lU ~l ' ' ' ; 3 

10 10 10 

Radius of gyration 

Fig. 3. Scaling of the number of broken bonds with the radius of the crack. 

In fig. (3) we show the typical behavior of the number of broken bonds. As an 
example we show data for a = 0.98 in which we averaged over four cracks. 
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Fig. 4. Dependence of fractal dimension on the breaking threshold. 

One finds that this number has a power law dependence on the radius of gyration 

N (X Rq' (11) 

and for this specific example we find a fractal dimension Df = 1.15. As indicated 
above this exponent varies with varying a. So, in fig. (4) we show all exponents Df 
for different a. 

Also in this plot the exponents Df{a) were obtained by averaging over four 
independent realizations. One obtains a linear dependence 

Df{a) = -(2.06 ± 0.05) • a + (3.19 ± 0.05). (12) 

Thus, for decreasing a one approaches a space filling structure which one reaches for 
a K, 0.58. However, it could be possible that for even larger clusters the asymptotic 
behavior changes and one crosses over into other exponents Df. 

6. Conclusions 

A'V'c! have described the implementation and results of a discrete time simulation for 
the growth of large cracks on a triangular network. We use a central force model 
and study the dynamical behavior of crack growth instead of studying the slow 
growth modes. Using a simplified dynamics, which anyway reproduces the essential 
features for the crack production, we arc able to grow cracks with more than 10000 
sites on a 1024 x 1024 lattice. We obtain fractal growth patterns with dimensions 
almost in the whole range between 1.1 and 2.0. The dimensionality of the cracks is 
mainly governed by the breaking threshold a and one finds a linear dependence of 
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